Method and apparatus for dynamic magnetic field control using multiple magnets

ABSTRACT

A method and apparatus for dynamic magnetic field control using multiple magnets. Control methods and system means are described that allow dynamically changing the magnetic field generated at a point in space by a multiplicity of magnets.

CROSS-REFERENCE TO RELATED APPLICATIONS

This application claims the benefit of U.S. Provisional PatentApplication Ser. No. 60/706,990 filed on Aug. 10, 2005, the entiredisclosure of which is incorporated herein by reference.

FIELD

This invention relates to the dynamic control of a magnetic fieldgenerated by two or more articulated magnets.

BACKGROUND

The robotic control and articulation of objects has been the subject ofmany investigations. The control of the direction and magnitude of amagnetic field created by a multiplicity of articulated source magnetsis a more complex topic that has been studied only more recently. Insuch a situation, the relatively simple geometric relationships definingtransformations between fixed and moving coordinate systems used inordinary robotics must be supplemented with knowledge of the magneticfields generated by the source magnets. These fields are typicallyrepresented by non-linear functions that most often are nottrigonometric in nature. While this problem has been investigated in afew specific cases, no general method has been described that allowsautomatic control of the magnets to generate a specific magnetic fieldat a given point in space. Embodiments of the present invention describea general solution to this problem.

One of the embodiments of the methods disclosed applies tooperator-directed control of the magnetic field direction and magnitudegenerated by two articulated magnets in an operating region. In aspecific application, the magnets considered are permanent magnetsgenerally facing one another and the operating region is centered inbetween the two magnets. The operator-directed control is executed bymeans of electromechanical positioners that move the two magnetsappropriately. The six degrees of freedom (three for each magnet) alloweach magnet to translate along the line between the two magnet centers,and each magnet to rotate in two angles θ and φ along respective axes.Other similar embodiments that could be used to control a multiplicityof magnets to generate a specific magnetic field (direction andmagnitude) at a specific point in space are also included in the presentinvention.

In magnetic navigation applications, it is desirable to change theorientation (and possibly magnitude) of the magnetic field at theoperating point to orient the tip of a catheter or guide wire withrespect to the field. The catheter or guide wire is then guided throughthe subject's vasculature while being simultaneously pushed at theproximal insertion point. It is further desirable to provide a “properturn” from the initial to the final field vector at the given operatingpoint. A proper turn is a rotation of the magnetic field direction inthe plane formed by these two vectors. The concept of a proper turn andits execution is described in U.S. Pat. No. 6,702,804, issued Mar. 9,2004, from Application Ser. No. 09/678,640 filed Oct. 3, 2000, forMethod for Safely and Efficiently Navigating Magnetic Devices in theBody, incorporated herein by reference.

SUMMARY

The procedure to execute a proper turn in a preferred embodiment of thepresent invention entails setting a series of steps for the turn, i.e.finding intermediate field vectors (directions and magnitudes) whichcouple the two known (initial and target) field vectors at the operatingpoint through the turn.

For each intermediate magnetic field vector the problem of finding thearticulations or values of the parameters describing each of the degreesof freedom for the magnets necessary to generate the given field is theinverse problem of finding the field position, direction, and magnitudeof magnets having specific articulations.

For each intermediate magnetic field vector at the operating point, andin one embodiment of the present invention, the inverse problem referredto above is solved in three steps. In the first step, the magnitude ofthe desired field at the operating point is apportioned to the twocontributing fields of the two articulated magnets. In the second step,an approximate representation for each magnetic field generated by eachmagnet is used to ensure quick convergence of the algorithm to anapproximate solution. In the third step, optimization methods areapplied to a more accurate and mathematically complex fieldrepresentation to obtain a refined estimate of the magnets controlparameters and a more accurate field solution.

In the first step, the magnetic field apportionment algorithm operatesunder various constraints related to the nature of magnet design andfield distributions to optimize field uniformity near the operatingpoint.

In the second step, and in one embodiment of the present invention, anested polynomial representation for each of the magnet fields enablesan efficient solution to the inverse problem. These representations arebased on calibration data acquired prior to any navigation procedure.Optimization operates separately for each magnet in consideration of theapportioned fields.

In the third step, and given that the output of the second step is inthe neighborhood of the optimal solution, the inverse problem isformulated as a simultaneous optimization in all degrees of freedom, andthe approximate polynomial field representation of step two is replacedby a more accurate field representation. One embodiment of the presentinvention uses a spherical harmonics expansion field representation.

In accordance with one aspect of the present invention, a multiplicityof N magnets (N>=2) is controlled to achieve a target magnetic field ata point in space (the operating point). The control algorithm proceedsin three steps. In the first step, the target field is apportionedbetween the N magnets. In the second step, the inverse problem offinding the magnets articulations that provide the apportioned field atthe operating point is solved independently for each of the N magnetsusing a first level of magnetic field representation. In the third step,the problem of finding the magnet articulations that generate the targetfield at the operating point is formulated as a simultaneous constrainedoptimization problem for all N magnets. The third step uses the outputof step two as a starting point in the neighborhood of the optimalsolution and also uses a more accurate, second level of magnetic fieldrepresentation for optimization.

The magnets preferably comprise one or more three types: permanent,permanent focusing, or electromagnetic, and are preferably permanentmagnets.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 is a schematic representation of one of the magnets andassociated motion axes, including rotation angle θ and pivoting angle φ.

FIG. 2 presents a typical magnet field distribution and associatedmagnitude contours in the plane Y^(M)=0 (in the coordinate systemassociated with the magnet).

FIG. 3 shows the parameterization (ρ, α, β) for the orientation andmagnitude of the projection of the magnetic field B onto a given (X^(M),Z^(M)) plane.

FIG. 4 shows maximum Z values achievable as a function of the pivotingangle φ. This represents the spatial constraint imposed by the fixedmagnet covers.

FIG. 5 presents the navigation algorithm flowchart.

FIG. 6 illustrates the magnetic field vector progression from itsinitial state B₀ to its target state B^(T) at the operating pointthrough a succession of intermediate field values.

FIG. 7 shows the time progression of one of the state vector elements,Z, through a series of intermediate discrete values.

FIG. 8 presents the three-step inversion algorithm flowchart.

FIG. 9 illustrates for the case N=2 the constraints on each of the twofield magnitudes and a method to select apportioned values for the twomagnets. FIG. 9-A presents the case when the two magnitudes can be setequal; FIG. 9-B and FIG. 9-C illustrate alternative magnitudeapportionment choices when selecting equal magnitudes for the two fieldsdoes not lead to a satisfactory solution due to minimum and maximummagnitude constraints on each of the fields.

FIGS. 10-A and 10-B present respectively angle and magnitude errordistributions achieved by step 2 of the inverse algorithm, using anested polynomial representation for each of the fields and independentoptimization for each magnet.

FIGS. 11-A and 11-B present respectively corresponding angle andmagnitude error distributions achieved by the optimized inversealgorithm of step 3, using a local spherical harmonics functionexpansion and simultaneous optimization in all variables.

Corresponding reference numerals indicate corresponding parts throughoutthe several views of the drawings.

DETAILED DESCRIPTION

In one embodiment of the present invention, control of two permanentmagnets enables dynamic generation of a specific magnetic field at anoperating point. The operating point may be chosen anywhere within asphere of specific radius centered mid way between the two magnets intheir calibration frame positions (at θ_(i)=φ_(i)=0; i=1,2). To relatethe magnetic fields generated by each of the two magnets to that presentat the operating point, three Cartesian coordinates systems are defined:one associated with each of the two magnets, and a third associated withthe subject (world coordinate system). A given vector, such as thatrepresenting the desired magnetic field at the operating point withinthe subject, can be represented in each of these coordinate systemsthrough the usual coordinate transformations. The generation of adynamic magnetic field at the operating point entails motion of the twomagnets with respect to each other and also with respect to the subject.The position of the operating point in each magnet coordinate systemwill change, as well as the direction of the desired and intermediatefield vectors. Accordingly, a series of intermediate magnetic fieldvectors that lie in a plane in the subject coordinate system will not ingeneral lie in a plane of either of the magnet coordinate systems. Thelocation and orientation of the two magnets at a given time is calledthe state of the system at that time. The system state evolves in timeas necessary to carry out a proper turn.

The Cartesian referential frame (X^(M), Y^(M), Z^(M)) associated witheach magnet is described in FIG. 1. FIG. 2 illustrates the fielddistribution in the central (X^(M), Z^(M)) magnet plane. By design thisis a plane of even field symmetry with respect to Y^(M). Proceeding inthe Y^(M) direction, succeeding (X^(M), Z^(M)) planes have similar butgradually changing field patterns.

Returning to FIG. 1, and considering one magnet among the multiplicityof magnets under control, two mechanical rotations and a translationdefined with respect to the magnet coordinates allow for changing thefield distribution in the world (subject) coordinates. These motions aredefined as follows:

-   -   “Pivoting” in the text refers to a rotation by an angle φ, 1,        about the Y^(M) axis 2. This motion alters the contributions of        the X and Z components of this magnet to the field at the        operating point.    -   “Rotation” in the text refers to a rotation by an angle θ, 3,        about an axis that coincides with the Z^(M) axis 4 when the        pivot angle is zero; this motion controls the central symmetry        plane of the magnet, and is used to align this particular magnet        central plane with the target field B^(T) at the operating        point. The effect of this motion is to provide a Y component of        the field of this magnet at the operating point.    -   “Translation” in the text refers to translating the magnet by a        distance Z along an axis that coincides with the Z^(M) axis 4        when the pivot angle is zero; this motion controls the field        magnitude.        Functionally, rotation is used to orient the magnet central        plane(X^(M), Z^(M)) such that it contains the target field        B^(T). Pivoting adjusts the direction of the magnet axes (X^(M),        Z^(M)) so that axis −X^(M), 6, is aligned with B^(T) and the        Z^(M) component of B^(T) in the magnet central plane is reduced        to zero; if that component (the Z^(M) component of B^(T)) is        small to start with then the need for pivoting is reduced, as is        desirable from a mechanical stand point.

The following description of the preferred embodiment of the presentinvention is divided in three parts: A) Field representation; B)Navigation method; and C) Inversion method. In part A), various methodsof describing the magnetic fields at several levels of accuracy andcomplexity are described. Part B) further formulates the problem ofmagnetic navigation as that of fitting intermediate field vectors thatdefine a proper turn while meeting a number of design constraints. InPart C), for each of the intermediate magnetic field vectors, theinverse problem of finding the magnet articulations that generate such amagnetic field at the operating point is solved in three steps. Thefirst step apportions the total magnetic field magnitude to each of themagnets. Step two finds an approximate solution by solving the inverseproblem for each of the magnets independently, given the first level ofmagnetic field representation introduced above and the apportionedfields of step one. Step three finds an improved solution bysimultaneously solving the problem for all magnets and using a moreaccurate field representation.

Part A): Field Representation

The magnet design is meant to generate a specific magnetic field in athree dimensional volume. Independently of the magnet designspecifications, the generated fields can be represented to any level ofaccuracy by a finite element model or by three dimensional look-uptables. Such representations are inefficient in that they requirestorage of a large amount of data (in function of the desired accuracy)and lead to the use of computationally intensive three-dimensionalinterpolation methods. In one embodiment of this invention, the magneticfields requirements are specific to magnetic navigation of catheters andguide wires in the body of a subject. In this embodiment, the actualfields generated by magnets designed per the correspondingspecifications can be approximated to a very good accuracy by a basisexpansion on spherical harmonic functions using a finite number ofterms. In this representation, and referring to FIG. 1, the even and oddsymmetries of specific field components with respect to the magnet'sY^(M) and Z^(M) axes 2 and 4 are taken into account. Commonly, theexpansion coefficients are determined through a calibration procedureusing the field values on an imaginary sphere containing the magnet. Inparticular, such a representation can satisfy the Laplace equation thatapplies to stationary magnetic fields, and can also be used to locallyrepresent the field in a specific volume (such as within a sphereenclosing the operating point) to a very high degree of accuracy.

However, such an expansion on spherical harmonics functions may not bethe most efficient means to represent the field distribution within agiven accuracy and a limited number of coefficients. Further, such arepresentation might not be optimal in terms of performing the inversecalculations that are necessary to determine the system state that willgenerate a specific field at the operating point. In one embodiment ofthe present invention, these two limitations are addressed by definingtwo levels of field representation. The first level of accuracy usesnested polynomials to provide a simplified field representation thatallows efficient inversion algorithm convergence for each magnetindependently to a state vector in the neighborhood of the optimum. Inthe second level of field representation accuracy, the coefficients ofthe spherical harmonics expansion are calculated by performing a leastsquares minimization between the representation and the actual field ina specified volume of interest. The volume of interest can be defined bya sphere that will encompass the operating point under all practicalmagnetic navigation situations. Such a local representation allows botha minimized error within the volume of interest for a given number ofexpansion terms as well as improved convergence of the inversionalgorithm.

For a given magnet, on the central, symmetry plane (X^(M), Z^(M)), theY^(M) component of the magnetic field is small (nominally zero bydesign). Because the field direction and magnitude vary slowly withdistance e from the central plane, and because in a typical navigationconfiguration the Y^(M) field components of the two facing magnets tendto combine destructively (the respective y axes being coaxial but inopposite directions), in a first approximation it is natural ignore the) field components. Any actual state vector optimization result obtainedunder this assumption for a target field within the navigation volumewill include a small e field resultant. Under this assumption, the fieldfor each magnet at any point within a volume of interest can then berepresented by its orientation in the associated (X^(M), Y^(M)) planeand magnitude. Referring to FIG. 3, and introducing cylindricalcoordinates (ρ, α, γ), ρ radial distance 8, α radial angle 10 withrespect to axis Z^(M) 12, and γ elevation, the field at any point in themagnet navigation volume is uniquely determined by the knowledge of theangle β=β(ρ,α,γ) 14 and field magnitude |B|=|B(ρ,α,γ)|. A convenientrepresentation for β and |B| uses nested polynomials as described below.

This section details the field magnitude representation. The Bio-Savartlaw calls for the field magnitude created by an elemental electricalcircuit to vary according to the inverse of the square of the distance.It suggests representing the inverse of a macroscopic magnet fieldmagnitude by a polynomial function of distance ρ of limited degree. Asinversion of the problem is a necessary step of the method, it isdesirable to limit the polynomial order to two; such a representation isalso sufficiently accurate for the object of the method. Accordingly,the field magnitude on any point of an iso-|B| curve 16 of FIG. 3 ismodeled by:$\frac{1}{B} = {{a(\alpha)} + {{b(\alpha)}\rho} + {{c(\alpha)}{\rho^{2}.}}}$In a given (X^(M), Z^(M)) plane, the field Z-component presents an oddsymmetry with respect to axis Z^(M), and the variations along α for thecoefficients a, b, and c can be efficiently represented by afourth-degree polynomial in α:a(α)=a ₀ +a ₁ α+a ₂α² +. . . +a ₄α⁴b(α)=b ₀ +b ₁ α+b ₂α² +. . . +b ₄α⁴c(α)=c ₀ +c ₁ α+c ₂α² +. . . +c ₄α⁴As indicated above, the field variations along Y^(M) are relativelyslow, and accordingly each of the 15 coefficients a₀, . . . ,a₄,b₀, . .. ,b₄,c₀, . . . ,c₄ may in turn be represented by a fourth-degreepolynomial function of Y^(M). The resulting set of 75 coefficientsaccurately describes the field magnitude at any point in a volume thatencompasses the navigation sphere.

This section details the field angle representation. Given the magnituderepresentation |B|=|B(ρ,α,γ)| which for fixed α and γ is a one-to-onefunction of the distance ρ, the field angle β^(M) 14 can be representedby a function of |B| and α. That is, on a line 16 of constant |B|, thefield angle β^(M), defined as the angle between the projection of thefield vector onto the plane at Y^(M)=constant and the negative X axis ofthe (X^(M), Y^(M)) frame, is given by:β^(M)=ƒ(|B|, α).As will be further described below, the control algorithm considers anoperating point located on the Z^(M) axis. For a given apportionedtarget field B^(T) at the operating point, the projection of B^(T) onto(X^(M), Z^(M)) forms an angle β^(T), 0≦β^(T)≦π, with respect to the−X^(M) axis. The inverse algorithm then finds the angle α^(T) such thatβ^(M)=ƒ(|B|,α^(T))=β^(T)−α^(T). As FIG. 2 shows, β(0)≈0 and β growsmonotonically from β(−π/2)≈−π to β(π/2)≈π and therefore the functionβ(α) is monotonic and one-to-one over the interval [−π/2, π/2] and suchan angle can always be found. Consequently, after pivoting the magnet byan angle α=α^(T), the magnet field angle, β′=β^(M)+α^(T) (calculatedwith respect to the fixed world −X axis) will be equal to the targetangle β^(T). To simplify the calculations during navigation andapplication of the algorithm, it is therefore desirable to represent thefunction α(β) inverse of the one-to-one function β^(M)=ƒ(|B|, α). Withα(0)=0 a cubic polynomial fit can be written as:α=d(|B|)β+e(|B|),β³.The coefficients d and e can be modeled as quadratic functions of thefield magnitude:d(|B|)=d ₀ +d ₁ |B|+d ₂ |B| ²e(|B|)=e ₀ +e ₁ |B|+e ₂ |B| ²Each coefficient d₀, d₁, d₂, e₀, e₁, e₂, is then modeled as a cubicfunction of Y^(M). Accordingly these 24 coefficients provide an accuraterepresentation of the angle α(β).

These nested polynomial field representations for field magnitude andangle yield maximum errors of:

-   -   ε(|B|)≦1.5 mTesla, and    -   ε(β)≦0.7 degrees    -   for the range of variables relevant for navigation in the        subject volume.        Part B): Navigation Method

The navigation algorithm solves the problem of finding a time sequenceof source magnets states such that the final sequence state provides anaccurate and precise estimate of the target magnetic field at theoperating point. As described above, the state of the system iscompletely determined given the two magnet Z_(i) translationcoordinates, rotation angles θ_(i), and pivoting angles φ_(i) (i=1,2).The problem can be more formally stated as:

-   -   Given:    -   an operating point at r, and    -   x₀ (an initial system state) corresponding to an initial field        B₀ at the operating point, and    -   B^(T), the target field at the operating point,    -   Find x(t) the state evolution to obtain the target field at the        operating point.        It is desired that the magnetic field changes in a plane in the        subject frame as the system evolves from the initial state        x₀=x(t₀) to the final state x(t_(f)). In a normal navigation        procedure it is also desirable to keep the field magnitude        essentially constant; in a well designed magnet system the        relative field magnitude variation between initial and final        states is not expected to be larger than a few percents.        Constraints are imposed on the velocity and acceleration to        which each magnet may be subjected. Z translation constraints        are imposed by the fixed magnet covers and other design        considerations; also the magnet pivoting angle is limited to a        specific range. Mathematically, the procedure to find the        desired progression of state x(t) is then:        Minimize the departure of B(t) from the plane formed by B₀,        B^(T);        Minimize the difference between the field magnitude |B(t)| and a        functional of the initial and final magnitudes |B₀| and |B^(T)|        (such as average of |B₀| and |B^(T)| or a polynomial fit),        Subject to:    -   B_(f)=B(t_(f))=B^(T)    -   Z<Z_(n) (Z_(n) is the closest position to the system cover)    -   −40°<φ<40°    -   a_(min)<a<a_(max), for all variables Z, θ, φ (acceleration        conditions)    -   v_(min)<v<v_(max), for all variables Z, θ, φ (intermediate        velocity conditions)    -   v₀=v_(f)=0    -   x₀=x(t₀)    -   x_(f)=x(t_(f)).

Additional requirements are imposed on the control method. It is oftendesired in a magnetic medical procedure that the resultant fieldmagnitude at the operating point be kept essentially constant; as foreach magnet the contributed field magnitude varies with distance fromthe magnet center, this is accomplished by translating the magnets toalleviate the field magnitude changes that are unavoidably associatedwith magnet reorientation (rotation and pivoting). In the preferredembodiment of a magnetic navigation system fixed and close-fittingcovers enclose each of the magnets to ensure safety and subject access;the associated limits on magnet translation are magnetorientation-dependent, as illustrated in FIG. 4. The maximum Z valuealong axis 4 is shown as a function of the pivoting angle φ, 1. Incertain orientations the magnet must be pulled back a small distance toavoid collision with the cover. A number of other conditions are imposedto make the motions and articulations more efficient and safe; asexamples, it is desirable to limit both the accelerations to which themagnets are subject as well as the resulting velocities; mechanicalconsiderations indicate that it is also desirable to reduce the pivotingmotion (angle φ in FIG. 1). Generally it is favorable to apportion thefield magnitudes equally between the two magnets as this fielddistribution locally minimized field gradients. Yet other conditions andconstraints will be shown to be useful in the control method andalgorithms that enable performing a proper turn.

A large number of inversion and optimization approaches are described inthe literature. Linear problems formulated as least-squares estimationcan be solved using matrix approaches. Such problems are written as:Min∥Ax−b∥²Degenerate inversion problems (that is, problems with an infinite numberof solutions) can be regularized by including a weighted penalty term inthe cost function; the additive combination of a degenerate and anon-degenerate quadratic being non-degenerate.Non-linear constrained optimization problems have been investigated in anumber of settings. For a problem formulated as:Min ƒ(x)s.t. g _(i)(x)≦0,iεIa general formalism leading to practical solutions uses a Lagrangianfunction:${L\left( {x,\lambda} \right)} = {{f(x)} + {\sum\limits_{i \in I}{\lambda_{i}{g_{i}(x)}}}}$Necessary conditions for local optimum have been given by Kuhn andTucker. These conditions have also been extended to non-linear problemswith both equality and non-equality constraints. It is also clear fromthe formalism above that a linear problem, such as that of fitting apolynomial to a series of data point, under specific equality andinequality constraints, can be represented and solved through theLagrangian formalism.

The navigation algorithm proceeds in the following five steps, describedin more details below, and illustrated in the flowchart 40 of FIG. 5.Given B₀, B^(T), 44:

1. Define intermediate target field values at predetermined angles, 46;

2. For each intermediate field value, calculate by inversion theassociated state vectors (these intermediate states will then befitted), 48;

3. Assess the intermediate interval times based on velocity andacceleration constraints, and determine an initial total time, 50;

4. For a given total time 52, iterate on the following till all theinequality constraints are met or an upper bound on the number ofiteration is reached:

5. For the four end constraints (initial and end state and velocityconstraints) and additional weighted cost functions (as necessary torepresent specific navigation constraints), solve the optimizationproblem using a Lagrangian approach 54. If, 56, one or more of theinequality constraint(s) is not met, increase the weight on thecorresponding constraint, 58, and iterate. If the maximum number ofiteration has been met, 60, increase the total time allotted to themotion and iterate. The final state vector sequence 64 is given at thecompletion of the iteration.

For each magnet and for each of the magnet coordinates Z, θ, φ, therequirements described above translate into four equality constraints(two for the initial and end velocities, and two for the beginning andend state conditions). The inequality requirements translate intoadditional constraints (four to six depending on the variable). It ispossible to design exact methods that will achieve a proper turn byensuring that intermediate field values lie exactly in the plane formedby the initial and end field vectors. However such methods arecomputationally expensive and are subject to oscillations (as attempt atpolynomial interpolation will demonstrate). A more efficient methodconsists of specifying intermediate field values (themselves prescribinga proper turn), calculating by inversion the corresponding state values,and fitting these values by a functional. The inversion procedure givesthe intermediate state values to be fitted, rather than interpolated, bya polynomial. As there are four equality constraints, and additionalinequality constraints, a polynomial of at least degree six isprescribed. At least seven fit nodes or intermediate magnetic fieldvalues are used in this embodiment to keep the field close to that of aproper turn.

FIG. 6 shows schematically for a navigation turn the intermediate nodesto be fitted during navigation from an initial field B₀ 70 to a finalfield B^(T) 74 at the operating point 76. The intermediate field valuesB_(i), 72, are provided in 10-degrees increments.

The algorithm begins by solving the inverse problem for eachintermediate field value B_(i). The inverse algorithm, described belowin Part C), provides the intermediate state vector values to beapproximated during navigation through the fitting procedure. Next, thetotal time T_(f) for the motion is estimated as follows. For each axisand for each segment joining two adjacent state values the absolutecoordinate increment is calculated, and the total distance betweeninitial and final state vector coordinates is then estimated as the sumof the absolute (linear) increments over all segments. Assumingparabolic velocity profiles and using velocity and accelerationconstraints, and based on the total distances just calculated, thetravel time needed for each axis is calculated. The maximum of all thetotal axes motion times (plus a small extra amount) gives the totalinitial motion time T_(f). To find the intermediate times ti, the timesneeded for each axis to go from one intermediate position to the nextwith its maximum velocity are determined. The time needed by the slowestaxis in each interval is retained, and all the intermediate times arescaled so that their sum equals the total motion time T_(f).

An illustration of the resulting steps in magnet translation coordinateZ_(i) for this procedure is shown in FIG. 7. The time series 80 of Zvalues from Z₀ 82 to Z_(f) 84 for one magnet shows some fluctuation thatare to be expected. The polynomial fitting of these points is part ofthe polynomial method associated with this one of the six magnetcoordinates. Comparable figures would show the progression of steps in θand φ for each of the two magnets. A sixth degree polynomial suitablyfits the data. For Z(t) the time progression is modeled as (for thespecific basis function retained):Z(t)=y ₀ +y ₁ t+. . . +y ₆ t ⁶and similar equations are written for the variables θ(t), φ(t), for bothmagnets. It is noted that the t_(i) values are not necessarily equallyspaced. Solution for the polynomial coefficient can be carried out bymeans of the associated Vandermonde matrix (p=6):1 1·1t₁ t₂·t_(n)····t₁ ^(P) t₂ ^(P)·t_(n) ^(P)Other basis function choices would lead to other matrix forms. It isdesirable to scale the time values to a limited interval to increasenumerical stability. This can be done, for example, by normalizing thetime sample values by use of the series distribution mean m, andstandard deviation σ_(t): t_(i)=(t_(i)−m_(t))/σ_(t).

The optimization in terms of the y coefficients that fit the Z(t_(i))values must be made subject to the four end point constraint in x and v.A suitable cost function will need four components with relative weightschosen as follows:

1) Fit errors-largest weight;

2) Acceleration root-mean-square-second largest weight;

3) Z and φ root-mean-square distance from center of range ofmotion-third largest weight;

${\sum\limits_{i}y_{i}^{2}} - {{smallest}\quad{weight}}$(Vandermonde inversion regularization).The problem is then solved analytically for the y_(i) coefficients byminimizing the cost functions subject to the four end point constraintsv₀=v_(f)=0, x₀=x(t₀), x_(f)=x(t_(f)). This is solved by using aLagrangian function L, setting ∂L/∂x and ∂L/∂λ to 0 and finding thecorresponding Lagrange multipliers λ_(j).This procedure provides an initial set of fit coefficient values y₀, y₁,. . . , y₆ for which the end point constraints have been satisfied. Thenext algorithm step tests whether the inequality constraints are alsosatisfied. If not, the algorithm increases the weight on thecorresponding constraint(s). This procedure is iterated until allconstraints are satisfied. Finally, if any velocity or accelerationconstraint(s) is still violated, the time points are scaled, therebyincreasing the total motion time, and the procedure is iterated.

The same method can be used if B is increased or reduced in magnitudeonly, but in that case, the intermediate target field vectors can beinterpolated directly with no need for a polynomial fit.

Part C): Inverse Method

The inverse algorithm proceeds in three steps, further described below,and illustrated in the flowchart 90 of FIG. 8. In the first step 94, thetarget magnetic field magnitude b=|B^(T)| is apportioned between themultiple magnets. In one embodiment, the field magnitude b isapportioned between two magnets: b=b₁+b₂. Two methods are described forthis apportionment, a geometric method, and an optimization method thatalso enables consideration of specific constraints via the inclusion ofcorresponding cost function terms. In the second step 96, first theinitial rotation and pivoting are performed to i) align the desiredB_(i) ^(T) vector in the X^(M), Z^(M) plane and ii) to put the operatingpoint on the Z^(M) axis in order to determine the magnitude bounds foreach magnet. Then, using the nested polynomial field magnitude and anglerepresentations described above in Part A), the pivoting angle andcorresponding Z translation (if necessary) are found that solve thefield problem for each magnet independently. In the third step 98, themore accurate field representation using spherical harmonicrepresentation is used to solve the problem formulated as a simultaneousoptimization in all magnet coordinates and obtain an improved solution100 that more closely match the target magnetic field in angle andmagnitude.

Part C), Step One: Field Apportionment

This section details the geometric field magnitude allocation method.This method provides an apportionment of b=|B^(T)| to b₁, b₂ (andtherefore, as will be seen below, of B^(T) to B₁ and B₂) by consideringthe magnitude bounds for each magnet. During navigation, and to maintainfield strength, the two fields contribute constructively (additively,i.e., with both fields in the direction of the desired vector). In theprocedure known as field reduction (not performed during navigation of amedical device in a subject), the fields contribute destructively (i.e.,with the two fields in opposite directions).

For any particular target field, it is necessary to find whether theproposed apportioned fields are within the upper and lower bounds ofeach magnet. Different situations are illustrated in FIGS. 9-A through9-C, all showing plots of the dashed line b₁=b₂ 112 and the solid lineb=b₁+b₂ 114 with respect to the axes b₁ 116 and b₂ 118. The lower andupper bounds of the respective magnet fields B₁ ^(M) and B₂ ^(M) areindicated by t₁, 120, u₁, 122, l₂, 124, and u₂, 126. Intersection of thedashed line b₁=b₂ and the solid line b=b₁+b₂ provide appropriatecategorization. If b/2 is within the bounds of B₁ ^(M) and B₂ ^(M), thenb₁=b₂=b/2. This case is illustrated in 9-A, 130, where the intersectionof the two lines occurs within the rectangle 132 representing the magnetbounds. As described previously, this is a most favorable condition, asit minimizes the magnetic field gradient. For some field directions andsome locations within the navigation volume, these bounds will not bemet. FIGS. 9-B and 9-C show two cases for which the line crossings falloutside the rectangle formed by the magnets bounds. In such cases theapportionment is determined as follows: Take the intersection of therectangular region of valid magnitudes with the line b=b₁+b₂. If thisintersection is empty, then select the vertex of the rectangle closestto the line. If the intersection is non-empty, then find on thisintersection the point closest to the crossing of the two lines b=b₁+b₂and b₁=b₂. In FIG. 9-B, 140, the lines intersect at point A=(A₁, A₂),142, and the closest bound is at A′=(A′₁, A′₂), 144. Thus b₁=l₁=A′₁,120, and b₂=A′₂, 148, are the retained apportionments. In the case ofFIG. 9-C, 150, the lines cross at C, 152, and the closest valid point isC′=(C′₁,C′₂), 154,, so the corresponding apportionments areb₁=u₁=C′₁,122, and b₂=C′₂, 158.

This section details the penalized optimization field magnitudeallocation method. The problem is formulated as a constrainedoptimization problem. There are six variables, consisting of the threefield components for each of the two magnets. There are three equalityconstraints from the equation B₁+B₂=B. In addition there are fourinequality constraints from the bounds l₁≦b₁≦u₁ and l₂≦b₂≦u₂. Thecomponents of the cost functions are chosen as follows:

1) |b₁−b₂|², to reduce the field gradient;

2) [B_(1,Z)/b₁]²+[B_(2,Z)/b₂]², to reduce the fields components alongthe coordinate axis perpendicular to the magnet surface; this constraintminimizes the amount of pivoting required;

3) b₁ ²+b₂ ², to force the fields to add constructively.

The respective terms of the cost functions are combined by weightsdetermined by experience. The problem is solved in closed-form usingLagrange multipliers. A solution can always be found, considering thesmall number of constraint surfaces. There is a sufficiently smallnumber of Kuhn-Tucker points (stationary points of the Lagrangian) thateach satisfies first-order necessary conditions for a (local) minimum,that they can all be checked and compared in turn to determine theglobal minimum. A general expression for the Lagrangian is:L(x,λ,μ)=ƒ(x)+λ^(T) g(x)+μ^(T) h(x),where g(x) represents equality conditions and h(x) represents inequalityconditions.

In this section, the problem of apportioning a target magnetic fieldB^(T) to N magnets is considered:$B^{T} = {\sum\limits_{i = 1}^{i = N}B_{i}}$The geometric apportionment method for N=2 extends to the general casein a straightforward manner. As before, select the field direction foreach magnet the same as B^(T). In order to apportion the magnitudes,define the hyperplane${H = \left\{ {{\sum\limits_{i = 1}^{i = N}b_{i}} = b} \right\}},$the line L={b_(i)=b_(j), i_(j), =1, . . . , N}, the rectangularpolyhedron R={l_(i)≦b_(i)≦u_(i)}, and the point d=H∩L, L={b_(i)=b/N i=1,. . . , N}. If the intersection of H and R is empty, meaning that thereis no valid selection of individual field magnitudes that add up to thetarget magnitude, then select the vertex of R closest to H in Euclideannorm as the apportionment point. Since none of the edges of R can beparallel to H, this vertex is unique. If the intersection of H and R isnonempty, then it will define a polytope P of dimension N−1. Theapportionment point is selected as the point in P closest to d inEuclidean norm. Since P is a closed convex set, this point is unique.Computationally this point can be found by solving the following linearleast-squares problem:min{∥y−d∥², y εP}and setting b_(i)=y_(i). Since this is a quadratic program with linearequality and inequality constraints, it can be solved efficiently usingthe algorithm described in below.

Alternatively, the optimization field magnitude allocation method ofabove can be generalized for N magnets as follows: There are 3Nvariables, consisting of the three field components of the N magnets.There are three equality constraints from the equation${\sum\limits_{i = 1}^{i = N}B_{i}} = {B^{T}.}$In addition there are 2N inequality constraints from the boundsl_(i)≦b_(i)≦u_(i) i=1, . . . , N. The components of the cost functionsare chosen as follows:${\left. 1 \right)\quad{\sum\limits_{i = 1}^{i = 1}{{b_{i} - {b/N}}}^{2}}},$to reduce the field gradient;${\left. 2 \right)\quad{\sum\limits_{i = 1}^{i = N}\left( {B_{z,i}/{B_{i}}} \right)^{2}}},$to reduce the fields components along the world coordinate axisperpendicular to the magnet surface; this constraint minimizes theamount of pivoting required;${\left. 3 \right)\quad{\sum\limits_{i = 1}^{i = N}b_{i}^{2}}},$to force the fields to add constructively. This can be solved using aniterative algorithm for constrained nonlinear programming.Part C), Step 2: Independent Magnet Control using Polynomial FieldRepresentations

In the previous step, the field apportionment of b to b₁ and b₂ and thusof B^(T) to B₁ and B₂ was found. The problem solved in this section isto that of finding the state sub-vector x₁ that yield B₁ at theoperating point. This step is carried out independently for each magnetusing a first order magnet field representation.

First, the rotation angle θ₁ is found such that B₁ is contained in an(X^(M), Z^(M)) magnet M₁ plane. This is achieved by projecting B₁ ontothe (X^(M), Y^(M)) plane. The angle −θ₁ is then the angle between −X^(M)and the projection of B₁.

Second, the pivoting angle φ₁ is found that aligns the magnet Z^(M) axiswith the line from the magnet coordinate center to the operating point,i.e. brings the operating point to the (Y^(M), Z^(M)) plane. At thisstage, the field magnitude bounds are determined.

The vector B₁ now lies in an (X^(M), Z^(M)) magnet plane; B₁ forms aknown angle β₁ with the −X^(M) axis.

The next step is then to determine the combination of variables (φ₁, Z₁)that will provide a field matching B₁ at the operating point. This canbe achieved in two ways:

-   -   i) By performing bisection searches using a look-up table field        representation of the field (knowing that both the field        magnitude and the field angle are monotonic functions of the        distance p and of the pivoting angle φ). The use of such look-up        tables requires time-consuming three-dimensional interpolations;        or:    -   ii) In a preferred embodiment, by using the nested-polynomial        field representation previously introduced. In this        approximation, the magnet field component along Y^(M) is        ignored. First the pivoting angle α₁ is estimated by plugging β₁        into the functional α(β) and then the distance ρ is found by        inverting the quadratic form for 1/b. In this way, the target        position that gives the desired field B₁ is found in the magnet        frame using cylindrical coordinates. Finally, the magnet        translation Z (function of ρ) and pivoting angle φ (function of        α) are found from the system geometry and corresponding motions        are implemented to bring the actual target point to the desired        target position.        Part C), Step 3: Constrained Simultaneous Optimization using a        more Accurate Field Representation

In step 3, the problem of finding a state vector x that yields a totalcontributed field B(x) equal to the target field B^(T) at the operatingpoint is formulated as a penalized optimization problem by consideringthe multiple magnets simultaneously and by using a more accuratespherical harmonics field representation. This step will in generalimprove upon the solution found in step 2 and provide an optimizedsolution that more closely match the target magnetic field in angle andmagnitude.

As a result of step 2, a state vector estimate x is available; thatestimate was derived from independent inversion for each of the magnetsusing nested-polynomial field representations. The algorithm nowproceeds to find a state vector x (for all magnets) that minimizes theEuclidian distance between the target and the estimated fields|B^(T)−B(x)|², using a more accurate field representation. There arepotentially an infinite number of solutions. The Euclidian distance costfunction may have a long shallow valley, for which convergence isdifficult in at least one dimension. The following constraintsregularize the problem and allow quick convergence to an optimum. First,a component is added to the cost function to minimize the distance tothe solution of step 2 (with a small weight). Adding this term cups theend of the valleys, making the cost function locally convex andimproving the likelihood of finding a unique solution. Another,independent non-linear constraint is to avoid interference of thestationary covers with the moving magnets. This might require that asthe magnet pivots, the magnet be pulled back in Z. This constraint isshown in FIG. 4. The pivoting angle is also limited to the range:−40°<φ<40°, and Z>−4″. The rotation angle is not constrained. Thisoptimization problem is a nonlinear least-squares problem with upper andlower bounds on the variables and the nonlinear constraints due to thecovers. The problem is formulated as a penalized optimization problem,and the total cost function can be written as the summation of squaresof nonlinear functions ƒ _(i)(x). This problem is solved by thefollowing iterative procedure.

At each step, the cover-related nonlinear constraints are locallyapproximated with linear chords of the Z(φ) curves as exemplified inFIG. 4. The chord is chosen around the current point, in a range thatdepends on whether the curve is locally concave or convex. If duringalgorithm iterations the cover constraint is slightly violated due tothe chord approximation, the Z variable is adjusted to meet theconstraint. Each term ƒ_(i) squared in the cost function is approximatedby its linear term, thus resulting in a quadratic sub-optimizationproblem for this step. A search direction (also known as theGauss-Newton direction) d is selected from the linear manifold of activeconstraints. Along this direction, a line search using quadraticpolynomial fit approximation to the non-linear total cost function isperformed that yields a local minimum and the next state vector iterate.At the end of the step, a new constraint may become active, an formerlyactive constraint may have become inactive, or a global minimum point isreached, which is then the solution of the main optimization problem.

This algorithm has been tested by computer simulations for 30,000 pointsin the subject operation volume, in which the target field magnitude was0.1 Tesla. Step 2 yields angle errors of less than 1.7 degrees in 95% ofthe samples, with a maximum of 11 degrees. In comparison, step 3 reducesthe angle errors to less than 0.73 degrees in 95% of the samples, with amaximum of 2.9 degrees. For magnitude errors, step 2 yields less than5.3 mTesla in 95% of the samples with a maximum of 9.4 mTesla. Incomparions, step 3 reduces the magnitude errors to less than 4.3 mTeslain 95% of the samples, with a maximum of 8.4 mTesla. Histograms of theseresults are shown in FIGS. 10-A and 10-B for step 2 and FIGS. 11-A and11-B for step 3. On average, steps 1 and 2 together take 3 ms, and step3 takes 140 ms on a personal computer.

The extension of the methods of steps 2 and 3, described in sections[0050]-[0055] to more than two magnets (N≧3) is straightforward. Sincestep 2 is carried out for each magnet independently, it applies withoutmodification. In step 3, the same cost function and constraints can bewritten in a space of dimension 3N, and the iterative solution proceduredescribed above will find the solution, although it will take morecomputation time to perform.

Quasi-Continuous Navigation

In one embodiment of the algorithm, a target field B^(T) is generated bythe operator using an input device such as a joystick. At the same time,the operating point may remain fixed, or it may also be changed by theoperator using another input device, such as a pointing device, or itmay be automatically calculated by the system using a mechanical modelof the interventional device being navigated. A new target field or anew operating point may be requested while the magnets are moving, alongthe trajectory calculated by the algorithm, corresponding to thepreviously requested field. In this case, the algorithm is modified asfollows: The magnets are decelerated to stop at an intermediateposition; during the deceleration, the path calculation is repeated,this time starting from the intermediate position and ending at theposition corresponding to the newly requested target field; and themotion is started again. This may be repeated as often as necessary,whenever a new field request is generated by the operator, implementinga version of quasi-continuous navigation.

Haptic Feedback

During quasi-continuous navigation, it is useful to provide hapticfeedback to the operator by applying a resistive force to the joystick.The force will indicate the subjective difficulty of the requested fieldchange. The difficulty may be represented by making the resistive forceproportional to either one of the following alternatives, or acombination of these:

-   -   i. The weighted sum of the difference between the current and        the final joint positions, i.e. the distance traveled by each        joint.    -   ii. The estimated time required to reach the final joint state.        The time may be calculated as the maximum of the required times        for each joint using the distance traveled and a parabolic        velocity profile; or it may be calculated using the path        generation algorithm in Part B.    -   iii. The vector difference between the current and final field        vectors. This would provide a directional resistive force.        More Degrees of Freedom—Secondary Mechanism

More degrees of freedom may be added, by using a secondary articulationmechanism, in order to move each magnet together with its primaryarticulation mechanism, so as to obtain a better imaging angle or forpatient access etc. The secondary mechanism may be floor- orceiling-mounted, and may be constructed independent of the primarymechanism. During the motion of the secondary mechanism, therelationship between each magnet and its primary articulation mechanismremains the same, so that all the calculations described before for thethree axes of motion, Z, θ and φ apply without modification. However,the relationship between each magnet and the patient changes. Therefore,it is desirable to coordinate the motions of the primary and thesecondary mechanisms so that the generated field B remains the same inthe patient. This may be accomplished by modifying the algorithm asfollows: First intermediate points are selected along the motiontrajectory for the axes of the secondary mechanism. At each intermediatetime point, the transformation between the patient and each magnet frameis recalculated. Next, at each intermediate point, Part C Steps 1 and 2are repeated for each magnet to find the corresponding joint positionsof the primary mechanism to generate the target field. At the finaltime, Part C Step 3 is performed combining all magnets. Finally, Part Bis performed separately for each magnet to determine the motiontrajectory.

In one embodiment, the secondary mechanism performs a tilting motionthat rotates each magnet and its primary mechanism around the imagingisocenter without changing the relationship between the magnet frames.

Feedback from a Magnetic Sensor

If it is desired to improve the accuracy of the generated magnetic fieldbeyond the accuracy of the numerical representation in Part A, then onemay use a high precision magnetic field sensor to measure the actualfield vector B^(m) at the operating point. First, the navigationalgorithm described in Parts B and C is performed to bring the magnetsto a penultimate joint position x, with the corresponding measured fieldB^(m)(x). Let e(x)=B^(T)−B^(m)(x) be the vector error between the targetand measured fields. Denote the derivative of the field with respect tothe joint state vector x by A(x)=∂B(x)/∂x. Corresponding to a smallincrement v in the joint position vector, the linear approximation ofthe field around the penultimate position x is B(x+v)≈B(x)+A(x) v. Theincrement is calculated as the minimizer of the cost function∥A v−e∥²+v^(T)R vwhere R is a diagonal matrix with positive entries chosen by experienceas weights to penalize large changes in the corresponding jointvariable. Since the matrix R is positive definite, the unique incrementv is found by solving the linear system(A ^(T) A+R)v=A ^(T) e=A ^(T)(B ^(T) −B ^(m))which constitutes a feedback law using the magnetic sensor. Then, theincrement in the Z positions of the magnets are modified if necessary tosatisfy the cover constraint in Part B. The time required to performthis increment is determined as the shortest time that satisfies thevelocity and acceleration constraints in Part B.

If, after such an increment, the magnitude of error in the field vectoris larger than a specified value, then the increment procedure isrepeated starting from the modified joint vector x+v, until the error issmall enough or a given maximum number of iterations is reached

The advantages of the above described embodiment and improvements shouldbe readily apparent to one skilled in the art. Accordingly, it is notintended that the invention be limited by the particular embodiment orform described above, but by the appended claims.

1. A method of determining the positions and orientations of at leasttwo magnets to apply a magnetic field at an operating point of specifiedstrength and direction, the method comprising: (1) apportioning thecontribution to the desired field of each of the at least two magnets;(2) individually determining for each of the at least two magnets, aposition and orientation to achieve the apportioned contribution using arepresentation of the magnetic field of each the magnet; (3)simultaneously determining a position and orientation of each of the atleast two magnets two magnets to achieve the desired field, using theindividually determined positions and orientations and a representationof the magnetic field of each magnet.
 2. The method of claim 1 whereinthe representations of the magnetic field used in step (3) are differentfrom the representations of the magnetic field used in step (2).
 3. Themethod of claim 1 wherein step (2) comprises using a nested polynomialrepresentation of the field of each magnet.
 4. The method of claim 1wherein step (3) comprises using a spherical harmonics expansion fieldrepresentation of the field of each magnet.
 5. The method of claim 1wherein step (3) is performed as a simultaneous optimization of alldegrees of freedom.
 6. A method of generating a specific target magneticfield at an operating point in space, the method comprising the stepsof: (a) representing magnet and application related constraints; (b)defining a series of intermediate field vectors at the operating point;(c) using an inversion method to find a series of method controlparameters for a multiplicity of magnets associated with the series offield vectors of step (b); (d) finding a time sequence of controlparameters for the multiplicity of magnets of step (c) subject to theconstraints of step (a) to approximate the series of control parametersof step (c); and (e) performing the time sequence of step (d), wherebythe target magnetic field at the operating point is achieved at a givenlevel of accuracy and within a specific time interval.
 7. The method ofclaim 6, wherein the series of intermediate field vectors at theoperating point lie in the plane formed by the initial and the targetmagnetic field vectors at the operating point. 8-12. (canceled)
 13. Themethod of claim 6, wherein the time sequence is obtained by fitting theseries of magnets control parameters of step (c) by representing thedistance between the fit functions and the series as one of theLagrangian terms.
 14. A method of representing the magnetic fieldgenerated by a magnet, the method comprising the steps of: (a) providinga least a first level of magnetic field representation such that themagnet control parameters that will provide a specific magnetic fieldgenerated by the magnet at an operating point in space within a firstlevel of accuracy can be found efficiently; and (b) providing at least asecond level of magnetic field representation such that the magnetcontrol parameters that will provide a specific magnetic field generatedby a multiplicity of magnets at an operating point in space at a second,higher, level of accuracy can be found efficiently; and (c) designingthe magnetic field representations of steps (a) and (b) such that thecontrol parameters found by using the representations provided in step(a) for the multiplicity of magnets under consideration in step (b) willprovide an initial point in the neighborhood of the final solution forthe procedure of finding the control parameters using the representationprovided in step (b); whereby the control parameters for a multiplicityof magnets that will generate a specific magnetic field at an operatingpoint in space within a second, higher level of accuracy subject todesign and application related constraints can be found efficiently. 15.The method of claim 14, wherein the step (a) of providing at least afirst level of magnetic field representation comprises the use of nestedpolynomials.
 16. The method of claim 15, wherein the nested polynomialsprovide a one-to-one correspondence between a couple consisting of amagnetic field magnitude and magnetic field orientation within a magnetplane of constant coordinate Y^(M) and a triplet of magnet articulationsconsisting of a rotation angle θ, a pivoting angle φ, and a translationdistance Z
 17. The method of claim 14, wherein the step (b) of providinga second level of magnetic field representation comprises using aspherical harmonics function expansion, and calculating the expansioncoefficients for optimal accuracy given a number of expansion terms byminimizing the representation error in a local volume; said local volumealways containing the operating point.
 18. The method of claim 1 whereinthe apportioning comprises distributing a target magnetic field at anoperating point in space among the at least two magnets subject todesign and application related constraints, whereby the sum of themagnetic fields contributed by the at least two magnets at the operatingpoint is equal to the target magnetic field at the operating point. 19.The method of claim 18, further constraining the fields generated byeach controlled magnet to be aligned with the target magnetic field andsuch that the target field magnitude is distributed in equal parts toall controlled magnets.
 20. The method of claim 19, where adjustmentsare made to the field apportionments to account for each magnet lowerand upper field magnitudes at the operating point. 21-24. (canceled) 25.The method of claim 6 wherein the inversion method comprises the stepsof: (a) using a magnetic field representation; (b) representing magnetdesign and application related constraints; (c) finding the magnetcontrol parameters using the representation of step (a) and under theconstraints of step (b) that will provide a magnetic field at anoperating point in space, whereby a target magnetic field at theoperating point in space is obtained at a given level of accuracy. 26.The inversion method of claim 25 using a nested polynomial magneticfield representation and further comprising the steps of: (a) pivotingthe magnet by an angle p such that the operating point lies on the Z^(M)magnet axis and with a plane of constant Y^(M); (b) calculating therotation angle θ such that the apportioned target magnetic field lies ina plane of constant Y^(M) in the magnet under consideration; (c)calculating the target magnetic field angle β with respect to the magnet−X^(M) axis; (d) calculating the cylindrical coordinate angle α of apoint within the magnet such that the associated magnetic fieldgenerated by the magnet lies at an angle α+β with respect to magnet axis−X^(M); (e) calculating the cylindrical coordinate ρ of a point withinthe magnet at an angle α+β such that the magnet field magnitude at thatpoint equals that of the apportioned field; (f) calculating from planargeometry the magnet translation Z and pivoting φ that will bring thetarget point to the desired position within the magnet for generation ofthe apportioned field at the operating point; and (g) carrying out therotation of angle θ, pivoting of angle φ, and translation along Z^(M) tobring the target point to the desired point within the magnet. 27-28.(canceled)
 29. The inversion method of claim 25 using a sphericalharmonics expansion representation for the magnetic field of each magnetand further comprising the step of formulating the problem as anon-linear constrained optimization method.
 30. The inverse method ofclaim 29, further comprising the steps of: (a) at each algorithmiteration, locally modeling non-linear cost functions via their linearapproximations and including the square of the linear terms in asub-problem optimization cost function; (b) adding a quadratic term tothe cost function that varies according to the square of the distancebetween the current iteration state vector x and the initialoptimization state vector x₀; (c) solving at each iteration step thelocal quadratic sub-problem using known techniques of constrainedoptimization; and (d) iterating on steps (a) through (c) above till asolution to the initial optimization problem has been reached.
 31. Themethod of claim 6 wherein the inversion method for the multiplicity ofmagnets comprises the steps of: (a) using a multiplicity of magneticfield representations for each the magnets; (b) representing magnetdesign and application related constraints; (c) using a first inversionmethod in a first stage to find the magnets control parameters using afirst level of magnetic field representation of step (a) and under theconstraints of step (b) that will provide a magnetic field at anoperating point in space at a first level of accuracy; (d) using asecond inversion method in a second stage, the magnet control parametersof step (c) as a starting point, and a second level of magnetic fieldrepresentation of step (a), finding under the constraints of step (b)the magnet control parameters that will provide a magnetic field at theoperating point in space at a second, higher, level accuracy, wherebythe target magnetic field at the operating point in space is achieved ata high level of accuracy and within a specified time interval by themultiplicity of magnets under control. 32-53. (canceled)